Attractiveness of periodic orbits in parametrically 
forced systems with time-increasing friction 



Michele Bartuccelli 1 , Jonathan Deane 1 , Guido Gentile 2 

1 Department of Mathematics, University of Surrey, Guildford, GU2 7XH, UK 

2 Dipartimento di Matematica, Universita di Roma Tre, Roma, 1-00146, Italy 



E-mail: m.bartuccelli@surrey.ac.uk, j.deane@surrey.ac.uk, gentile@mat.uniroma3.it 



Abstract 



We consider dissipative one-dimensional systems subject to a periodic force and study numer- 
ically how a time-varying friction affects the dynamics. As a model system, particularly suited 
for numerical analysis, we investigate the driven cubic oscillator in the presence of friction. We 
find that, if the damping coefficient increases in time up to a final constant value, then the basins 
of attraction of the leading resonances are larger than they would have been if the coefficient had 
been fixed at that value since the beginning. From a quantitative point of view, the scenario 
depends both on the final value and the growth rate of the damping coefficient. The relevance of 
the results for the spin-orbit model are discussed in some detail. 

1 Introduction 

Take a one-dimensional system driven by an external force, 



with F a smooth function 27r-periodic in time t; here and henceforth the dot denotes derivative with 
respect to time. Then add a friction term. Usually friction is modelled as a term proportional to the 
velocity, so that the equations of motion become 



with the proportionality constant 7 referred to as the damping coefficient. 

As a consequence of friction attractors appear [33] . If the system is a perturbation of an integrable 
one, there is strong evidence that all attractors are either equilibrium points (if any) or periodic orbits 
with periods T which are rational multiples of the forcing period 2tt [5j. If 2tt/T = p/q, with p, q E IN 
and relatively prime, one says that the periodic orbit is a p : q resonance. For each attractor one 
can study the corresponding basin of attraction, that is the set of initial data which approach the 
attractor as time goes to infinity. If all motions are bounded, one expects the union of all basins of 
attraction to fill the entire phase space, up to a set of zero measure. This appears to be confirmed 
by numerical simulations [121 HH1 [TJ [5] . 

Recently such a scenario has been numerically investigated in several models of physical interest, 
such as the dissipative standard map [191 HH1 HQ] , the pendulum with oscillating suspension point [7], 



x + F(x,t) = 





7 > 




the driven cubic oscillator [5] and the spin-orbit model [2j [13] . What emerges from the numerical 
simulations is that, for fixed damping, only a finite number of either point or periodic attractors 
is present and every initial datum in phase space is attracted by one of them, according to Palis' 
conjecture |39l 140] . However, which attractors are really present and the sizes of their basins depend 
on the value of the damping coefficient. If the latter is very small then many periodic attractors can 
coexist. This phenomenon is usually called multistability; see \19\ I20| . By taking larger values for 
the damping coefficient, many of the attracting periodic orbits disappear and, eventually, when the 
coefficient becomes very large, only a few, if any, still persist: for every other resonance there is a 
threshold value for the damping coefficient above which the corresponding attractor disappears. 

In this paper we aim to study what happens when the friction is not fixed but grows in time, more 
precisely when the damping coefficient is not a constant but a slowly increasing function of time. 
This is a very natural scenario: it is reasonable to suppose that in many physical contexts dissipation 
tends to increase to some asymptotic value. We aim to show (numerically) in such a setting, with 
the damping coefficient slowly increasing to a final value, that the relative areas of some basins of 
attraction become larger than they would be if the damping were fixed for all time at the final value. 
In other words, we claim that, in order to understand the dynamics of a forced system in the presence 
of damping, not only is the final value of the friction important, but also the time evolution of the 
damping itself plays a role. So, by looking in the present at a damped system which evolved from an 
original nearly conservative one, with the friction slowly increasing from virtually zero to the present 
value, it can happen that an attractor, which should have a small basin of attraction on the basis 
of the final value of the friction, is instead much larger than expected. Of course, as we shall see, 
several elements come into the picture, in particular the growth rate of the damping coefficient and 
the closeness between its threshold and final values. 

We shall investigate in detail a model system, the driven cubic oscillator in the presence of 
friction, which is particularly suited for numerical investigations because of its simplicity. We shall 
first study in Section [2] some properties in the case of constant friction, with some details worked 
out in Appendix |A| Then in Section [3] we shall see how the behaviour of the system is affected by 
the presence of a non-constant, in fact slowly increasing friction. In Section [4] we shall introduce 
another system of physical interest, the spin-orbit system, with some details deferred to Appendices 
[B] to [E] and we shall discuss how the results described in the previous sections may be relevant to 
the study of its behaviour. Further comments are deferred to Section [5] Finally in Section [6] we 
draw our conclusions and briefly discuss open problems. Some discussion of the codes we used for 
the numerical analysis is given in Appendix |F| 

2 The driven cubic oscillator with constant friction 

Let us consider the cubic oscillator, subject to periodic forcing and in the presence of friction, 

x + x 3 + ef(t)x 3 + 7± = 0, /(i) = cosi, (2.1) 

where x £ K, and e is a real parameter, called the perturbation parameter, that we shall suppose 
positive (for definiteness) . Of course one could consider more general expressions for / and the 



choice made here is for simplicity. The system (2.1) has been investigated in [3], with 7 a fixed 



positive constant. The constants e and 7 are two control parameters, measuring respectively the 
forcing and the dissipation of the system. 
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We shall look at ( |2.1[ ) as a non- autonomous first order differential equation, so that the phase 
space is K 2 . Note that (x,x) = (0,0) is an equilibrium point for all values of e and 7. Moreover, 
for e = 7 = the system is integrable and all motions are periodic. One can write the solutions 
explicitly in terms of elliptic integrals [25, 5j. For e ^ (hence e > 0) and 7 > fixed, a finite number 
of periodic orbits of the unperturbed system persist and, together with the equilibrium point, they 
attract every trajectory in phase space [5]. Such periodic orbits are called subharmonic solutions 
in the literature Each periodic orbit can be identified through the corresponding frequency or, 
better, the ratio u) := p/q between its frequency and the frequency of the forcing term. For each 
periodic orbit one can compute the corresponding threshold value j(u),s): if 7 > j(u),e) the orbit 
ceases to exist, while for 7 < 7(0;, e) the orbit is present, with a basin of attraction whose area 
depends on the actual value of 7. At fixed e, one has 



lim j(p/q, e) 

max{p,g}— >oo 



0. 



(2.2) 



Therefore, if we assume that all attractors different from the equilibrium point are periodic and no 
periodic attractors other than subharmonic solutions exist, then we find that at fixed e and 7 only 
a finite number of attractors exists. We note that the assumption above, even though we have no 
proof, is consistent with numerical findings [5j. 

The threshold value 7(0;, e) depends smoothly on e |14l |26| I22j : for all w G Q there exists 
n(ui) G IN such that the corresponding threshold value is of the form 7(0;, e) = Cq{uj, e) e n<yU) \ with 
the constant Cq(oj, e) nearly independent of e for e small; more precisely Cq{uj, e) tends to a constant 
Cq{oj) as e goes to zero, so that we can consider it a constant for e small enough. Resonances are 
classified as follows: we refer to resonances with frequency u such that n(uj) = 1 as primary, to 
resonances with frequency oj such that n{oS) = 2 as secondary, and so on |43} [5]. Of course such a 
classification makes sense only for e small enough. The primary resonances are the most important, 
in the sense that, at fixed small e, for 7 large enough, only primary resonances are present; moreover, 
by decreasing the value of 7, although non-primary resonances appear, they have a small basin of 
attraction with respect to those of the primary ones. The threshold values of the leading attractors 
(that is the attractors with largest threshold values), in terms of the constants Cq{oj), were computed 



analytically in [5] and are reproduced in Tables 2.1 and 2.2 In particular the periodic attractors with 
frequency 1/q, with q odd, appear in pairs [5j. The higher order corrections to Co(u,e) are explicitly 
computable; however we shall not need to do this here. Note that the classification of resonances 
and the corresponding threshold values strongly depend on the forcing: all values in this and next 
Section refer to f(t) = cost, as in (2.1). 



Table 2.1: Values of the constants Co(p/q) for p = 1 and q — 2, 4, 6, 8, 10 (leading primary resonances) for the 
cubic oscillator (2.1); the threshold values are of the form 7(0;, e) = Co(u>)e + 0(e 2 ). 



q 2 


4 


6 


8 


10 


C (l/q) 0.178442 


0.061574 


0.008980 


0.000920 


0.000078 



Consider the system (2.1 ) at fixed e. For 7 large enough, the only attr actor left is the equilibrium 
point; in that case all trajectories eventually go toward this point, which becomes a global attractor 
(see Appendix [A]) . If 7 is not too large — that is, according to Table 2.1, if 7 < Co(l/2)e, up to 



3 



Table 2.2: Values of the constants Co(p/q) for p = 1 and q = 1, 3, 5, 7, 9 (leading secondary resonances) for the 
cubic oscillator (2.1); the threshold values are of the form 7(0;, e) = Cq(ui)s 2 + 0(e 3 ). 



1 1 
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5 


7 


9 


C (l/q) 0.146322 


0.065001 


0.006488 


0.000177 


0.000002 



higher order corrections — then, besides the equilibrium point, there is a finite number of other 
attractors, which are periodic orbits. 

For e = 0.1, from Tables 2.1 and 2.2 one obtains the threshold values 7(1/2,0.1) ~ 0.018, 
7(1/4,0.1) « 0.0062, 7(1,0.1) ~ 0.0015, 7(1/6,0.1) ss 0.00090, 7(1/3,0.1) rj 0.00065, 7(1/8,0.1) ss 
0.000092, 7(1/5,0.1) ~ 0.000065, and so on. Take the initial data in a finite domain of phase space, 
say the square Q = [—1,1] x [—1,1]: the relative areas of the parts of the basins of attraction 



contained in Q for some values of 7 are given in Table 2.3 In principle the relative areas depend on 
the domain, but one expects that they do not change too much by changing the domain, provided 
the latter is not too small. Note that for e = 0.1 and 7 < 0.00005, other attractors than those listed 
in Table 2.3 appear (namely periodic orbits with frequencies 1/8, 1/5 and 3/4 for 7 = 0.00005; with 

0.00001; and with frequencies 1/10, 1/8, 1/7, 
0.000005), so explaining why the relative 



frequencies 1/8, 1/5, 3/10, 2/5, 5/12 and 3/4 for 7 
1/5, 3/14, 2/7, 3/10, 2/5, 5/12, 3/7, 2/3 and 3/4 for 7 
areas of the basins of attractions considered there do not sum up to 100%. Small discrepancies for 
the other values are simply due to round-off error (the error on the data is in the first decimal digit; 
see Appendix [F]) . 

Table 2.3: Relative areas A(uj,j), %, of the parts of the basins of attraction contained inside the square Q for 
e = 0.1 and some values of 7. The attractors are identified by the corresponding frequency (0 is the origin). 
The number of random initial conditions taken in Q is 1000 000 up to 7 = 0.0001, 500 000 for 7 = 0.00005, 
150 000 for 7 = 0.00001 and 50 000 for 7 = 0.000005. 
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7 = 0.020000 
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00.0 


00.0 


00.0 
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7 = 0.015000 
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7 = 0.010000 
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7 = 0.005000 
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7 = 0.000050 
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02.8 
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01.7 
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00.6 


7 = 0.000010 
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7 = 0.000005 
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If one plots the relative areas A(uj, 7) of the basins of attraction versus 7 one finds the situation 
depicted in Figure 2.1 Of course in general A(u,^f) depends also on e, i.e. A(u, 7) = A(u,^,e), 
although we are not making explicit such a dependence since e has been fixed at e = 0.1; same 
comment applies to the quantities A max (u) = ^4 max (a;,e) and A(oS) = A(uj,e) to be introduced. 
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Figure 2.1: Relative areas A(u>, 7) of the basins of attraction versus log 7 for the values of 7 listed in Table 
(right-hand figure) and a magnification for the periodic orbits with uj = 1/4, 1, 1/6, 1/3, 3/8 (left-hand figure). 



It can be seen that for any uj £ Q one has A(uj,j) = if 7 > 7(0;, e). By decreasing 7 below 
~y((jj,e), A(u, 7) increases up to a maximum value A m3iX (uj) which tends to stabilise. For very small 
values of 7 one observes a slight bending downward. It would be interesting to investigate very — in 
principle arbitrarily — small values of 7, but of course we have to cope with the technical limitations 
of computation: studying arbitrarily small friction would require running programs for arbitrarily 
large times and with arbitrarily high precision — see also comments in Appendix [FJ However, by 
looking at Figure [27T] and noting that numerical evidence suggests that all attractors different from 
the origin are periodic orbits, we make the following conjecture: as 7 goes to + , for all uj £ Q the 
relative area A(uj, 7) tends to a finite limit A{uj) such that 



5>m 



100%, 



(2.3) 



wed 



where uj = designates the origin. Of course when 7 = the area of each basin drops to zero, so 
that, accepting the conjecture above, all functions A{uj,~j) are discontinuous at 7 = 0. This is not 
surprising: a similar situation arises in the absence of forcing, where the only attractor is the origin, 
with a basin of attraction which passes abruptly from zero (7 = 0) to 100% (7 > 0). Moreover, 
analogously to [35], we would expect that the total number of periodic attractors N p grows as a 
power of 7 when 7 tends to 0. This means that if the limits A(uj) vanished at least one function 



A(lj,j) should be exponential in log 7, a behaviour which seems unlikely in the light of Figure 2.1 



3 The driven cubic oscillator with increasing friction 

Here we shall consider 7 = j(t) explicitly depending on time, that is 

x + x 3 + ef(t)x 3 + j(t)x = 0, f(t)=cost, 



(3.1) 



5 



For both concreteness and simplicity reasons, we shall consider a dissipation ^(t) linearly increasing 
in time up to some final value, i.e 



lit) 



7o ; 



70, 



< t < T , 
t>T , 



(3.2) 



where the parameters 70 and To are positive constants. However, the results we are going to describe 
should not depend too much on the exact form of the function 7(i), as long as it is a slowly increasing 
function; see Section [5] In (3.2) we shall take To = A/70, whose form is suggested by the fact that 
trajectories converge toward an attractor at a rate proportional to I/70 (see Appendix |A|) . 

Hence, consider the system with e = 0.1 again but now with 7(f) given by (3.2), with To = A/70 
and 70 = 0.015. Computing the corresponding relative areas A(u,^/o;A) = A(u), 70, 0.1; A) — that 
is A(u), 70, e; A) for e = 0.1 — for different values of A, we obtain the results in Table 3.1 and Figure 
|3.1| If A is very small, the damping coefficient reaches the asymptotic value 70 almost immediately, 
and we would expect to obtain the same scenario as in the previous case (7 constant): two attractors, 
corresponding to the origin and the 1:2 resonance, with basins whose relative areas are close to the 
values for A = 0, i.e. 91.1% and 8.9%, respectively. On the other hand, if A becomes larger, we find 
that the relative area of the basin of attraction of the origin decreases, whereas that of the basin of 
attraction of the 1:2 resonance increases. For A very large, these areas apparently tend to constant 



values of around 61% and 39%, respectively; see Table 3.1 



Table 3.1: Relative areas A(ui, 0.015; A) of the parts of the basins of attraction contained inside the square Q 
for e = 0.1 and j(t) given by (3.2) with 70 = 0.015 and T = A/70, for various values of A and w 
(oj = is the origin). In each case, 1 000 000 random initial conditions have been taken in Q. 
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08.9 


29.4 
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Analogous results are found, for instance, for 70 = 0.005; see Table 3.2 and Figure 3.2 One 
sees that for A = 20 the relative areas of the basins of attraction of the origin and of the 1:2 and 
1:4 resonances have already appreciably changed: they have become, respectively, 53.4%, 38.5% and 
8.1%. By further increasing A, once again a saturation phenomenon is observed and the relative 
areas settle about asymptotic values around 45%, 42% and 13% (for instance for A = 120 the areas 
are, respectively, 46.0%, 41.3% and 12.7%). Note that the value 70 = 0.005 is such that the threshold 
values 7(1/2,0.1) w 0.018 and 7(1/4,0.1) w 0.0062 of the persisting resonances are slightly above it 
(that is their ratios with 70 are of order 1). 

If one fixes the value 70 = 0.0005, then the 1:6, 1:1 and 1:3 resonances are also present. On 
the other hand the threshold values of the 1:2 and 1:4 resonances are appreciably larger then 70 
(that is 7(0;) S> 70 for ui = 1/2 and oj = 1/4), whereas the threshold values 7(1/6,0.1) w 0.00090, 
7(1,0.1) w 0.0015 and 7(1/3,0.1) ps 0.00065 of the 1:6, 1:1 and 1:3 resonances, respectively, are not 
too different from 70. If we again take j(t) as in (3.2), with To = A/70 and 70 = 0.0005, we have the 



results in Table 3.3 and Figure 3.3 



Therefore we have, from a qualitative point of view, the same scenario as in the case 70 = 0.005, 



6 



100.0 




200 



Figure 3.1: Relative measures of the basins of attraction versus A for 70 = 0.015. 



Table 3.2: Relative areas A(ui, 0.005; A) of the parts of the basins of attraction contained inside the square Q 
for e = 0.1 and j(t) given by (3.2) with 70 = 0.005 and To = A/70, for various values of A and u = 0, 1/2, 1/4 
(oj = is the origin). 500 000 random initial conditions have been taken in Q. 



A 


20 


40 


60 


80 


100 


120 


uj = 


64.8 


53.4 


49.4 


47.9 


46.8 


46.1 


46.0 


uj = 1/2 


31.8 


38.5 


40.0 


40.7 


41.2 


41.5 


41.3 


uj = 1/4 


03.4 


08.1 


10.6 


11.4 


12.0 


12.4 


12.7 



but with some relevant quantitative differences: the relative areas of the basins of the 1:2 and 1:4 
resonances are not too different in the two situations 7 constant and 7 increasing. 

We now give an argument to explain why the basins of attraction are different if 7 is not fixed ab 
initio to some value 70 but slowly tends to that value. According to Table 2.3, for smaller values of 7 
the basins of the periodic attractors are larger. For instance for 7 = 0.005 the 1:2 and 1:4 resonances 
have basins with relative areas 31.8% and 3.4%, respectively, while the basins of attraction of the 
same resonances for 7 = 0.0005 have relative areas 41.8% and 14.7%, respectively. Then, if we 
suppose that the friction is slowly increasing in time, when it passes, say, from 7 = 0.0005 to 0.005, 
on the one hand the size of the basin would decrease because of the larger value of 7, but on the other 
hand many trajectories have already nearly reached the basin and hence continue to be attracted 
toward that resonance. If friction increases slowly enough we can assume that it is quasi-static. 
Therefore, at every instant r, the basin of attraction of any resonance has the size corresponding to 
the value 7(7") at that instant, as can be deduced by interpolation from Table 2.3 (or Figure ??), 
while the rate of approach to the resonance can be roughly estimated as proportional to 1/7(7"); see 
Appendix [A} Therefore if A is large enough (that is if the growth of 7(f) is slow enough) one expects 
the trajectory to be captured by the resonance faster than how the basin of attraction is decreasing. 

By increasing the friction further, the basin of a resonance uj can become negligible, until the 
resonance itself disappears. If this does not happen, that is if 7(0;, 0.1) > 70, then there is a value 
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Figure 3.2: Relative areas of the basins of attraction versus A for 70 = 0.005. 



Table 3.3: Relative areas A(lj, 0.0005; A) of the parts of the basins of attraction contained inside Q for e = 0.1 
and 7(t) given by p^J with 7o = 0.0005 and T = A/70, for various values of A and uj = 0, 1/2, 1/4, 1,1/6, 1/3 
(uj = is the origin). 250 000 random initial conditions have been taken in Q. 
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of A above which the relative measure A(o;, 7 o; A) of the basin saturates to a value close to the 
maximum value A max (u) (possibly a bit smaller because of the slight bending downward observed in 
Figure 2.1). In particular, this explains the difference between Figures 3.2 and 3.3 For concreteness 
let us focus on the 1:2 resonance. With respect to the case 70 = 0.0005, according to Figure 3.1 the 
relative area A(l/2, 7) of the basin of attraction does not increase appreciably when taking smaller 
values 7 < 70: indeed ^4(1/2,0.0005) is already close to -A max (l/2). 

We conclude that the main effect of friction slowly growing to a final value 70, is that eventually 
every basin of attraction has essentially the same size that would appear for lower values of friction. 
So, if the basin of attraction of any p : q resonance is larger for values of friction lower than the final 
value, then, when the final value 70 is reached, one observes a basin of attraction with relative area 
larger than A(p/q, 70). If on the contrary it is more or less the same, then one observes essentially 
the same basin one would have by taking the friction fixed at that value since the beginning. In other 
words, if the friction increases in time, one can really have a larger basin of attraction only if the 
final value 70 is close enough to the threshold value 7(0;, 0.1). However, if 70 is too close to 7(0;, 0.1), 
for the phenomenon to really occur, the rate of growth has to be slow enough: the closer 70 is to 
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Figure 3.3: Relative areas of the basins of attraction versus A for 70 = 0.0005. 



7(0;, 0.1), the larger A to be chosen. For instance, for e = 0.1 and 70 = 0.0005, a glance at Table 3.3 



gives the following picture. The areas of the basins of attraction for the 1:2 and 1:4 resonances have 
small variations for different values of A, whereas the areas of the basins of attraction for 1:1, 1:3 
and 1:6 change in a more appreciable way when A becomes larger. Moreover, the threshold value 
7(1/3, 0.1) « 0.00065 is just above 70, so for the area of the corresponding basin of attraction to come 
close to the maximum possible value one needs large values of A; on the contrary 7(1, 0.1) ~ 0.00146 
is not too close to 70 and hence the area of the corresponding basin of attraction comes closer to the 
maximum possible value for smaller A ~ 10. 

We finish this section with a pair of figures showing the difference between the basins of attraction 
of the 1:2 resonance for 7 = 0.005, in the cases of constant and time-varying 7. According to Table 
changing A from zero — i.e. constant 7 — to 40 increases the area by about 8%, and Figure 
shows how this extra area is distributed (most of the points of the basin of attraction of the 



3^2 
3.4 



resonance for constant 7 still belong to the basin of attraction for varying 7) . 



4 The spin-orbit model 

The spin-orbit model describes the motion of an asymmetric ellipsoidal celestial body (satellite) 
which moves in a Keplerian elliptic orbit around a central body (primary) and rotates around an axis 
orthogonal to the orbit plane |23} 136] . If 6 denotes the angle between the longest axis of the satellite 
and the perihelion line, in the presence of tidal friction the model is described by the equation 

+ eG(0,t) + 7(0-l) = O, (4.1) 



where £,7 > and 9 £ T = H/2tt'Z, so that the phase space is T x R (note that (4.1) is of the 



form (1.2), with x = — t). Here e is a small parameter, related to the asymmetry of the equatorial 
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Figure 3.4: Basins of attraction for the 1:2 resonance for constant 7 (red) and time- varying 7 (black plus most 
of the red region). Initial conditions in the white region either go to the origin or to the 1:4 resonance. The 
right-hand figure shows a portion of the left-hand figure, magnified. 



moments of inertia of the satellite, and G(0,t) = dgg(6,t), where 



9(0, t) 



+ 
+ 
+ 
+ 



1 



1 

32 



123 



-e H e 

4 32 



768 

3 



e 5 )cos(26»-t) + 



489 rA , nn . (11 n 115 



e 4 ) cos(2#-4t) 



845 o 32525 

e H 

96 1536 

228347 



V 32 



7680 



e b cos(20 - 7t) + 



96 



e 4 cos(2<9 - 6t) 



e 5 cos(2(9 + i) 



(4.2) 



11 

1536' 



48 



e 4 ) cos(26> + 2i) + 



81 
2560 1 



e 5 cos(20 + 3i), 



with e being the eccentricity of the orbit; terms of order 0(e 6 ) have been neglected. In the celes- 



tial mechanics cases the model (4.1) may appear oversimplified and more realistic pictures could be 



devised [321 ESI IS] • Nevertheless, because of its simplicity, it is suitable also for analytical investiga- 
tions (as opposed to just numerical ones) on the relevance of friction in the early stages of evolution 
of celestial bodies and for the selection of structurally stable periodic motions. 

Values of e, e and 7 for some primary-satellite systems of the solar system are given in Table 
4.1 The values of e can be found in the literature [36] . while the derivation of e and 7 is discussed 
in Appendix |B| see also below. All satellites of the solar system are trapped in the 1:1 resonance 
(rotation period equal to the revolution period), with the remarkable exception of Mercury (which 
can be considered as a satellite of the Sun), which turns out to be in a 3:2 resonance. The spin-orbit 
model has been used since the seminal paper by Goldreich and Peale [23] in an effort to explain the 
anomalous behaviour of Mercury. The ultimate reason is speculated to be related to the large value 
of the eccentricity. However, even though higher than for the other primary-satellite systems, the 
probability of capture into the 3:2 resonance is still found to be rather low [231 03] • I n t ne following 
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part of this section we aim to investigate what happens if we take into account the fact that friction 
increased during the evolution history of the satellites. 



Table 4.1: Values of the constants e, e and 7 for some cases of physical interest for the spin-orbit model (4.1 ) 





Primary 


Satellite 


e 


e 


7 


E-M 


Earth 


Moon 


0.0549 


6.75xl0" 7 


3.75X10" 8 


S-M 


Sun 


Mercury 


0.2056 


8.11X10" 7 


3.24X10" 8 


J-G 


Jupiter 


Ganymede 


0.0013 


4.30X10" 4 


1.91X10" 5 


J-I 


Jupiter 


Io 


0.0041 


3.85X10" 3 


1.71X10" 4 


S-E 


Saturn 


Enceladus 


0.0047 


1.41X10" 2 


6.26xl0" 4 


S-D 


Saturn 


Dione 


0.0022 


3.85X10" 3 


1.71X10" 4 



Again one can compute the threshold values of the primary resonances, by writing 7(0;) = Cq{uS)s, 
up to higher order corrections. If one writes (|4.2|) as 



g(9, t) = ^a k cos(26> - kt), 



(4.3) 



one finds (see Appendix [Pj) 

Co(p/q) 



2q\ a 2p/q\ 

\p-q\ 



V 



' 2 2 2 2 



(4.4) 



while Co(l) = 00 (that is no threshold value exists for the 1:1 resonance) and Cq(uj) = for any 
other oj; other resonances may appear only at higher order in e. This leads to the values listed in 



Table 4.2, for the primary resonances of the systems considered in Table 4.1 Note that for 7 large 
enough, all attractors disappear except the 1:1 resonance, which becomes a global attractor. 



Table 4.2: Values of the constants Co(p/q) for some primary resonances of the the spin-orbit model (4.1 ); the 
threshold values are of the form 7(0;, e) = Co(u)e. Only positive ui have been explicitly considered. 





0J 


1/2 


3/2 


2 


5/2 


3 


7/2 


E-M 




5.488X10" 2 


3.818X10" 1 


2.545x10-2 


1.928X10" 3 


1.513X10" 4 


1.186X10- 5 


S-M 




2.045X10" 1 


1.308 


3.251X10" 1 


9.163x10-2 


2.976x10-2 


8.739X10" 3 


J-G 


Cb(w) 


1.300X10" 3 


9.100X10" 3 


1.436X10" 5 


2.578xl0~ 8 


4.757X10- 11 


8.832X10- 14 


J-I 




4.100xl0~ 3 


2.870X10- 2 


1.429X10- 4 


8.O88XIO- 7 


4.707X10" 9 


2.756X10- 11 


S-E 




4.700xl0~ 3 


3. 290X10" 2 


1.878X10- 4 


1.218X10- 6 


8.128X10" 9 


5.455X10- 11 


S-D 


C (u>) 


2.200X10" 3 


1.540x10-2 


4.114X10" 5 


1.259xlQ- 7 


3.902X1Q- 10 


1.226X1Q- 12 



It seems reasonable on physical grounds (see Appendix [E]) to assume that friction was increasing 
in the past up to the present-day value 70 = 7, with 7 as in Table |4.1| Application to the spin-orbit 
model for S-M would require numerics with very small values of e and 70: the discussion in Appendix 
[B] provides the values in Table [41] so that 70 ~ 0.05 e. However, the value usually taken for e in the 
literature is e ~ 10 -4 (see Appendix pL In both cases, 70 is far below the threshold value 7(0;, e), 
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especially for the most interesting resonance u = 3/2 [231I36], as 7(3/2, e) ~ 1.3e (see Table 4.2); we 
refer to Section [6] for further comments. 

As already noted in |13| . the small value of 70 represents a serious difficulty from a numerical 
point of view, because it requires very long integration (for a very large number of initial conditions 
- see comments in Appendix [F]) . Nevertheless, the discussion in Section [3] about the driven cubic 
oscillator allows us to draw the following conclusions about the spin-orbit model. 

The results in Table 2.3 suggest that the relative areas A(uj, r y,e) of the basins of attraction are 
almost constant for values of 7 much smaller than the threshold values; more precisely they assume 
values close to A max (uj, e). In the case of increasing friction, if the final value 70 is much smaller than 
the threshold value, the basin turns out to have more or less the same size close to ^4 max (u>,£) as it 
would have if 7 were set equal to 70 since the beginning. The same scenario is expected in the case 
of the spin-orbit problem. In particular a crucial aspect is understanding when the final value can be 
considered 'much smaller' than the threshold value or comparable to it. Again the analysis in Section 
[3] is useful: pragmatically, we shall define 7 much smaller than 7(0;, e) when A(oj, 7, e) is close to the 
maximum possible value ^4 max (a;,e). Therefore it becomes fundamental to check whether, for the 
current values of 7 and e as given in the literature, A(u, 7, e) is either close to A ma _ x (uj, e) or much 
smaller. 



1. If we assumed A (a;, 7, e) to be close to A ma _ x (u,e) the time-dependence of friction would not 
change the general picture as observed today. From this point of view, our results would be 
a bit disappointing: indeed the relative area of the basin of attraction of the 3:2 resonance, 
with the values of e and 7 usually taken in the literature, is found to be rather small for 
S-M |13j and including the time-dependence of friction in the analysis would not give larger 
estimates. On the other hand, also in this second case, our analysis would provide some more 
information: it would yield that the results available in the literature |15|. \W[ [T3] would remain 
correct even if time-dependent friction were included. In particular, to explain why Mercury 
has been captured into the 3:2 resonance, other mechanisms should be be invoked, such as the 
chaotic evolution of its orbit |15| . Of course the values of e and 7 used in the literature are only 
speculative: again our analysis suggests that the results would not change in a sensible way 
even by taking different values for one or both parameters — see also comments in Section [6| 

2. On the contrary if A(uj, 7, e) were much smaller than ^4 max (u;, e), taking into account the time- 
dependence of friction would imply a larger basin of attraction with respect to the case of 
constant friction. In this case the the exact values of the parameters e and 7 would play a 
fundamental role — again see also Section [6| 



5 Remarks and comments 



5.1 Different values of the perturbation parameter 



In Sections[2]and[3]we have fixed e = 0.1. However, by taking different values of e, the phenomenology 
does not change. For instance, for e = 0.5 and e = 0.01 we have the relative areas listed in Tables 



5.1 and 5.2 (for e = 0.5 only a few attractors are taken into account in Table 5.1) 
The general scenario is the same as in the case e 



to the fact that for smaller values of e (say e 



0.1, with obvious quantitative differences due 
0.01) only primary resonances are relevant unless 7 is 
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Table 5.1: Relative areas of the parts of the basins of attraction contained inside the square Q for e = 0.5. 
(ui = denotes the origin). 1 000 000 random initial conditions have been taken in Q. 



oj 


1/2 


1/4 la 


lb 1/6 


7 = 0.1000 


100.0 


00.0 


00.0 


00.0 


00.0 


00.0 


7 = 0.0750 


70.5 


29.5 


00.0 


00.0 


00.0 


00.0 


7 = 0.0500 


49.8 


50.2 


00.0 


00.0 


00.0 


00.0 


7 = 0.0250 


32.0 


56.0 


05.8 


03.1 


03.1 


00.0 


7 = 0.0050 


10.4 


48.8 


07.7 


10.5 


10.5 


00.0 


7 = 0.0025 


08.1 


36.3 


06.8 


11.5 


11.5 


00.8 


7 = 0.0010 


06.9 


37.5 


04.4 


11.4 


11.4 


01.6 



Table 5.2: Relative areas of the parts of the basins of attraction contained inside the square Q for e = 0.01. 
(oj = denotes the origin). 500 000 random initial conditions have been taken in Q. 



OJ 





1/2 


1/4 


la 


lb 


1/6 


7 = 0.0020 


100.0 


00.0 


00.0 


00.0 


00.0 


00.0 


7 = 0.0015 


98.1 


01.9 


00.0 


00.0 


00.0 


00.0 


7 = 0.0010 


94.0 


06.0 


00.0 


00.0 


00.0 


00.0 


7 = 0.0005 


88.3 


10.6 


01.1 


0.00 


0.00 


00.0 


7 = 0.0001 


78.2 


15.3 


06.5 


00.0 


00.0 


00.0 



very small, while for larger e (say e = 0.5) more and more resonances appear by taking smaller values 
of 7 (because powers of e are not much smaller than e itself). Indeed, if e = 0.5, even for 7 = 0.005 
we detect numerically more than 50 attractors and the classification of resonance oj according to the 
the value of n(oj) becomes meaningless. For instance, for e = 0.5, one has A(3/4, 0.0025, 0.5) ~ 11.4% 
and A(3/4, 0.001, 0.5) ~ 6.8%, to be compared with the values for uj = 1/4 in Table 5.1 Moreover for 
large values of e, say for e = 0.5, the bending of the curves ^(7, uj, e) in Figure ?? is more pronounced 
and the monotonic decrease observed for e = 0.1 when 7 tends to seems to be violated (compare 
the values A(l/2, 7, 0.5) for 7 = 0.0025 and 7 = 0.001 in Table 5.1 ); see also the comments in Section 

El 



5.2 Different functions y(t) 

As stated at the beginning of Section |3j the exact form of the function ^(t) should not be relevant. 
As a check we studied (3.1) with j(t) given by both (|3.2[) and 



7 (i) = 7o(l -e" i/To ), 



(5.1) 



where 70 and To are positive constants, by setting To = A/70 and changing A for fixed values of 70. 
The results show that the same behaviour is obtained in both cases. For instance, for 70 = 0.006, 



one has the results in Tables |5.3| and |5.4| see also Figure |5.1| 

In particular, in both cases the relative areas A(uj, 70, 0.1; A) start at the values corresponding to 
constant dissipation 7 = 70 for A = and then either decrease (for the origin) or increase (for the 
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Table 5.3: Relative areas A(oj, 0.006, 0.1; A) of the parts of the basins of attraction contained inside the 
square Q for e = 0.1 and j(t) given by (3.2) with 70 = 0.006 and Tq = A/70, for various values of A and 
oj = 0, 1/2, 1/4. (lu — denotes the origin). 500 000 random initial conditions have been taken in Q. 



A 


10 


20 


30 


40 


50 


oj = 


69.6 


64.3 


56.3 


53.1 


51.4 


50.1 


oj = 1/2 


29.9 


34.8 


37.2 


38.6 


39.2 


39.8 


oj = 1/4 


00.5 


00.9 


06.5 


08.3 


09.4 


10.1 



Table 5.4: Relative areas A(lu, 0.006, 0.1; A) of the parts of the basins of attraction contained inside the square 



Q for e = 0.1 and j(t) given by (5.1 1 with 70 = 0.006 and T = A/70, for various values of A and lu = 0, 1/2, 1/4 
(oj = denotes the origin). 500 000 random initial conditions have been taken in Q. 



A 





1 


2 


8 


13 


20 


30 


40 


to = 


69.6 


69.3 


67.0 


58.3 


55.8 


53.3 


51.3 


50.0 


oj = 1/2 


29.9 


30.1 


31.8 


36.1 


37.1 


38.2 


39.1 


39.7 


oj = 1/4 


00.5 


00.6 


01.2 


05.6 


07.2 


08.5 


09.6 


10.3 



periodic orbits), apparently to some asymptotic value close to ^4 max (u;, 0.1). For instance the relative 
areas corresponding to (3.2) with A = 30 are very close to those corresponding to (5.1) with A = 20 
(compare Table 5.3 with Table 5.4). The asymptotic values seem to be the same in both cases. 



6 Conclusions and open problems 

In this paper we have studied how the slow growth of friction may affect the asymptotic behaviour 
of dissipative dynamical systems. We have focused on a simple paradigmatic model, the periodically 
driven cubic oscillator, particularly suited for numerical investigations. Nevertheless we think that the 
results hold in the more general setting considered in Section [TJ The main result, discussed in Section 
[3j can be summarised as follows: on the one hand it is the final value of the damping coefficient that 
determines which attractors are present, but on the other hand the sizes of the corresponding basins 
of attraction strongly depend on the full evolution of the damping coefficient itself, in particular on 
its growth rate. Let "f(p/q, e) and 70 denote the threshold value of the p : q resonance and the final 
value of the damping coefficient, respectively. If 70 > "f(p/q,£) the attractor disappears. Otherwise 
the following possibilities arise: if 70 <C lip/q-,^) the area of the corresponding basin of attraction 
is more or less the same as in the case with constant damping coefficient equal to 70, whereas it is 
larger if 70 < ^(p/q, £)■ In the latter case, the closer 70 is to ^(p/q, e), the slower the growth rate 
of j(t) required for the maximum possible area of the basin of attraction to be attained. Moreover 
when 70 <C ^(p/q, e) the area can be even a little smaller than what it would be if the damping were 
fixed at 70 since t = 0, because other attractors may have acquired a larger basin of attraction while 
the damping coefficient increased. It is not possible for the relative area to be larger than the value 
Ama,x(p/q, e), which therefore represents an upper bound. 

Finally, let us mention a few open problems which would deserve further investigation, also in 
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Figure 5.1: Relative areas of the basins of attraction versus A for j(t) given by (3.2 ), with 70 — 0.006 (left-hand 
figure) and for j(t) given by (5.1|, with 70 = 0.006 (right-hand figure). 



relation with the spin-orbit problem. 

1. As far as a dynamical system can be considered as a perturbation of an integrable one, all 
attractors seem to be either equilibrium points or periodic orbits: at least, this is what emerges 
from numerical simulations. It would be interesting to have a proof of this behaviour, even for 



some simple model such as (2.1), in particular of the fact that neither strange attractors nor 



periodic solutions other than subharmonic ones appear. 

2. The analysis presented in Section [3] covers small values of 7, but still not as small as desirable. 
It would be worthwhile to study the behaviour of the curves A{oj^,e) in the limit of even 
smaller 7, for instance by implementing some numerical integrator which allows us to decrease 
the running time of the programs without losing accuracy in the results. 

3. In studying the spin-orbit model in Section [4] if one really wanted to consider the past history of 
the system, then not only 7 but also e should be taken to depend on time: this would introduce 
further difficulties and a more detailed model for the evolution of the satellite would be needed 
(see also comments at the end of Appendix [B|) . Of course, this would be by no means an easy 
task. The use of a model for the time-dependence of friction already raises several problems, 
as we highlight in Appendix |E} 

4. We have seen that, for the spin-orbit model, in the case of constant friction the exact value of the 
parameters e and 7 is fundamental. For instance if 7 = 71 such that A(o;,7i,e) « ^4 max (w,e), 
then we have a basin of attraction much larger than for 7 = 72, where 72 3> 71 is close to 
the threshold value 7(0;, e). On the contrary, our analysis shows that, when one takes into 
account that friction slowly increased during the solidification process of the satellite, then 
for both values 71 and 72 we expect more or less the same relative area close to A mg , x (u),e). 
This is useful information because it shows that exact values of the parameters e and 7 are 
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not essential in the case of time-dependent friction, as long as 7 is not much smaller than the 
threshold value (we mean 'much smaller' in the sense of Section [4]): of course the values of 7 
and e are essential if 7 is much smaller than j(uj, e) — see next item. 

5. As remarked in Section[4j to study the probability of capture of Mercury into the 3:2 resonance, 
it becomes crucial to determine the value ^4 max (3/2, e) and the current values of 7 and e to see 
if 7 turns out to be much smaller than 7(3/2, e), that is if one has ,4(3/2, 7) <C ,4 max (3/2, s). If 
this were the case, by using time-dependent friction, the relative area of the basin of attraction 
of the 3:2 resonance would be much larger than what was found for constant damping coefficient 
7 (estimated around 13% in p3]). By assuming the values of 7 given in the literature, for which 
7 ~ 1(T 8 and e ~ 1CT 4 , and using that 7(3/2, e) ~ 1.3 s, the relation ,4(3/2, 7) < ^ max (3/2,e) 
could be verified only by assuming for ,4(3/2,7, e) a very slow variation in 7 for e ~ 10 . 
This does not seem impossible: already for the driven cubic oscillator the profiles of the curves 



A(u, 7, e) seem to have a much smoother variation for smaller values of e (compare Table 2.3 for 



e = 0.1 with Table 5.1 for e = 0.5), so it could happen then by decreasing e further the curves 
A{uj, 7, s) could be nearly flat on much longer intervals: in other words, by fixing e ~ 10 -4 and 
decreasing 7, the curve A(u), r y,£) could have not yet reached its maximum value ,4 max (u;,£) 
at 7 ~ 10 -8 . As a consequence, the exact values of the parameters e and 7 and the profile of 
the curve ,4(3/2, e, 7) could be fundamental. In any case, it would be very interesting to study 
numerically the spin-orbit model for very small values of e and 7, in order to obtain the profiles 
of the curves -4(3/2,7, e) versus 7. 

Acknowledgments. We are grateful to Giovanni Gallavotti for very useful discussions and critical 
remarks, especially on the spin-orbit model and the formation and evolution of celestial bodies. 



A Some analytical results on system (2.1) 
A.l Global attraction to the origin for large 7 

For e small enough introduce the positive function F(t) such that F 2 (t) = 1 + ef(t) and rescale time 
through the Liouville transformation |35|, |6] 

r:= I dsF(s). (A.l) 
Jo 



Then we can rewrite (|2.1|) as 

' x' =y, 



y ' = - x3 -W)^ + F(t) 
where the prime denotes derivative with respect to r. Define 



V ( , F'(t)\ (A.2) 



2 4 

I( x ,x,t):=^ + ^J> ( A - 3 ) 
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which is an invariant for (2.1 ) with e = 7 = 0. More generally one has 



1 



7 + 



so that, if 



2(1 + £/(*» 



7> - min 2(i";; /( t))' (A - 4) 

by using Barbashin-Krasovsky-La Salle's theorem |31| . we find that the origin is an asymptotically 
stable equilibrium point and every initial datum is attracted to it as t — > 00. 



A. 2 Rate of convergence to attractors 



The periodic orbits for the system (2.1) appear in pairs of stable and unstable orbits: this is a 
consequence of Poincare-Birkhoff 's theorem [9]. Let us consider the primary resonances, so that we 
can set 7 = eC, with C a constant independent of e. 

We rewrite the equations of motion in action-angle variables (I, ip) [5] 

(if = (3/) 1 / 3 + e(3/) 1 / 3 /(t)cnV - eCcnipsiupdmp, 
[ I = e(3/) 4 / 3 /(t) cn 3 if sn p dn ip - e C(3 J) sn 2 </?dn (p, 

where cinp,snip,dnip are the cosine-amplitude, sine-amplitude, delta-amplitude functions, respec- 
tively, with elliptic modulus k = 1/V2 |25j . 

Let K (k) be the complete elliptic integral of the first type. For e = the periodic solution to (2.1 ) 
with frequency u = p/q is of the form x(t) =acn (a(t + to)), with 2ira = 4wK(l/2) and to suitably 
fixed [5]. In terms of action-angle variable this gives I = Iq '■= a 3 /3 and <p = (po{t) := a(t + to). 



Linearisation of (A. 5) around the periodic solution leads to the system 



^ - L(t) ( S *) , L(t) = Lo + eL.it) + O(e'), (A.6) 



where 



with 



L u (t) = f(t) d(t) -Ca" 1 ^), Li 2 (t) = -2a" 5 /i(t) + a- 2 f(t) a(t), 

L 2 i(t) = a 3 f{t)c{t)-Ca 2 d(t), L 22 (t) = 4a f(t) c(t) - 3C d(t), (A.* 



where we have defined 



a(t) := cn 4 <p (t), &(*) := cn(f {t) sncp (t) do.(p (t), 

c(t) := cn 3 ipo(t) sncpo(t) dn.<po(t), d(t) := sn 2 (po(t) cn 2 (po(t) (A. 9) 

and denoted by (</?i(i), -fi(i)) the first order of the periodic solution. 



17 



Let us denote by W(t) = Wo(t) +eW\{t) + 0(e 2 ) the Wrons kian matrix, that is the matrix whose 
columns are two independent solutions of the linearised system (A. 2), so that W(t) = L(t)W(t), with 
W(0) = 1. Then one has 



Wo(t) = exptLc 



1 



v° 1 

while W\(t) is obtained by solving the system W x = LqW x + L x Wo(t), i.e. 



W^t) = W (t) 



W x (0) + f 'driWoir))- 1 L x (t)W (t) 
J o 



(A.10) 



(A.11) 



where one has to take Wi(0) = in order to have W(0) = 1. 



A trivial computation shows that in (A. 10) one has 
W Q (t)- l L x {t)W Q (t) = 



L u (t) - a-HL 21 {t) L 12 (t) + a~ 2 t(L n {t) - L 22 (t) - a~ 2 L 21 (t)) 
L 21 (t) L 22 {t) + crH L 21 {t) 

Let T = 2nq be the period of the periodic solution. The Floquet multipliers around the periodic 
solution are the eigenvalues of the Wronskian matrix, computed at time T. Denote by x k (t) the kth 
primitive of any function x(t) with x k (0) = (so that x k (t) = x k ~ 1 (t), with x°(t) = x(t)). Then, by 
using that 

f T f T 

/ dtx(t) = x 1 (T), / dttx(t) = Tx l {T) - x 2 (T), 
Jo Jo 

/ dt t 2 x(t) = T 2 x x (T) - 2Tx 2 (T) + 2x 3 (T) , 
Jo 

we obtain that 



(A.12) 



_ (L\ X {T) + a- 2 L 2 21 (T) L\ 2 (T) + a^iTL^T) - L\ x + L 2 2 + a' 2 {TL 2 21 - 2Ll l {T)))\ 

For e = 0, the corresponding Floquet multipliers are equal to 1. To first order they are the roots A-t 
of the equation A 2 — 26qA + cq = 0, with 



b := 1 + - (Ll^T) + L\ 2 {T) + a^TLl^T)) , c := 1 + e (L^T) + L l 22 (T)) 



so that 



One has 



A± = 1 ± Jea-*TL\y{T) + - [L\ X {T) + L l 22 {T) + a^TL^T)] + o{e). 



L\ X {T) + L\ 2 {T) = [ dtf{t)(a(t) + 4ac(t))-C f dt (a _1 6(t) + 3d(t)) 
io Jo 

One immediately realises that the first integral vanishes and hence 



L\ X (T) + 4 2 (T) = -3CTfi, 



1 

T 



T 



dtd{t) > 0, 



(A.13) 
(A.14) 
(A.15) 

(A.16) 
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while, for the stable periodic solution, to is such that L\ X {T) < 0. Therefore the Floquet multipliers 
are of the form 

A± = 1 ± iXoVTe - (Ag + SCfj) Te + o(e), 
with Aq > 0. The corresponding Lyapunov exponents, defined as T _1 Re log A±, are given by 



—3fiCe = —3/^7- This shows that for primary resonances of the system (2.1), at least for initial 
conditions close enough to the attractors, convergence to the attractors has rate 1/7. In principle 
the analysis can be extended to any resonance, by writing 7 = Ce m and going up to order m, for 
suitable m depending on the resonance (m = n(oj) for the resonance with frequency to; see Section 
[2]) : the contributions to the Lyapunov exponents due to the Hamiltonian components of the vector 
fields cancel out and the leading part of the remaining part turns out proportional to —7. 



In the case of the linearly increasing friction (3.2) one expects that the Lyapunov exponent be still 



proportional to —7. If the friction increases very slowly, one may reason as it were nearly constant 
over long time intervals, that is time intervals covering several periods, by approximating 7(f) with 
a piecewise constant function. For each of such interval 7 can be considered as constant and one can 
reason as above. When passing from an interval to another, the value of the initial phase to of the 
attractor slightly changes. The Lyapunov exponent is then expected to behave proportionally to 



1 r f T f T ° A 

lim- / dt 7 (t), / dt 7 (t)= / dt7(t)+7o(r-r ) = - + ( 70 r- A), 

t^oo T J Jo Jo 2 



where we have used that To = A/70. Therefore again the rate of exponential convergence to the 
attractor is proportional to I/70 and after the time To the distance to the attractor has already 
decreased by a factor exp(— coA) = exp(— C070T)), for some positive constant Co, and hence like in 
the case of constant friction, possibly with a different constant cq. 



B Parameters £ and 7 for the spin-orbit model 

The spin-orbit model has been extensively used in the literature to study the behaviour of regular 
satellites |23[ [36] — it does not apply to irregular satellites, which are very distant from the planet 
and follow an inclined, highly eccentric and often retrograde orbit. The equations of motion are given 

by 

9 + eG{0,t) = O, (B.l) 



with G as in (4.2). Here time has been rescaled t — > wyt, where ujt is the mean angular velocity of 



the satellite along its elliptic orbit (cf. Table B.l), so that the orbital period ('year') of the satellite 
becomes 1. Then the 1:1 resonance is 9 ~ 1. 

In a system satellite-primary there can be several types of friction: for instance the friction 
between the satellite layers of different composition, say one liquid and one solid (core-mantle friction), 
or the friction due to the tides (tidal friction). One can expect that such phenomena produce a friction 
to be minimised in a 1:1 resonance. There could be also other sources of friction which we do not 
consider, especially those which could modify the revolution motion of the satellite, because we are 
implicitly using that it occurs on a fixed orbit. The dissipation term due to tidal torques is of the 
form [33 [231 SH M 

- 7 (fi(e)0-JV(e)), (B.2) 
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Table B.l: Values of lot (angular velocity), M (satellite mass), Mq (primary mass), R (satellite radius) and p 
(mean distance between satellite and primary) for the systems considered in Section HI CGS units are used. 





Primary 


Satellite 




M 


M | R 


P 


E-M 


Earth 


Moon 


2.66xl0" 6 


7.35X10 25 


5.97X10 27 


1.74x10 s 


3.84X10 10 


S-M 


Sun 


Mercury 


8.27X10" 7 


3.30X10 26 


1.99x10 s3 


2.44x10 s 


5.79X10 12 


J-G 


Jupiter 


Ganymede 


1.02X10" 5 


1.48X10 26 


1.90X10 30 


2.63x10 s 


1.07X10 11 


J-I 


Jupiter 


lo 


4.11X1CT 5 


8.93X10 25 


1.90X10 30 


1.82x10 s 


4.22X10 10 


S-E 


Saturn 


Enceladus 


5.31X1CT 5 


1.08X10 23 


5.68X10 29 


2.52X10 7 


2.38X10 10 


S-D 


Saturn 


Dione 


2.66X1CT 5 


1.09X10 24 


5.68X10 29 


5.62X10 7 


3.77xl0 10 



where 0(e) and iV(e) are two constants depending on e. Since both O(e) 
1 + 0(e 2 ), for small values of e we can approximate (B.l) as in (4.1). 
A comparison with the literature [23, 15J gives 



1 + 0(e 2 ) and N(e) 



3/„ 



3/i 
2R' 



7 



£Q\p) V M ) ' 



(B.3) 



where I x ,I y ,I z are the moments of inertia of the satellite, h is the maximal equatorial deformation 
(tide excursion), R and M are the mean radius and the mass of the satellite, Mq is the mass of 
the primary, p is the mean distance between satellite and primary and &2, £, Q are constants, known 
respectively as the potential Love number, the structure constant and the quality factor. For instance, 
in the case of the Moon one has k 2 ~ 0.02, £ « 0.4 and Q w 30 [281 EE], which g ives 3 V£<2 ~ 0.005 
and hence 7 ~ 3.75 x 10 -8 (approximately 3.15 x 10 -6 years -1 ). In the case of S-M, the constants are 
usually (somehow arbitrarily) set equal to the values k 2 ~ 0.4, £ w 0.3333 and Q m 50 [25 } [2^1 15"4" 1 [2T] . 
which gives 3k2/£Q ~ 0.072. The corresponding value of the damping coefficient is 7 ~ 3.24 x 10 -8 , 
a value very close to Mercury's. Expressed in years -1 this becomes approximately 8.46 x 10 -7 (the 
revolution period of Mercury is 27r/o;r ~ 7.60 x 10 6 s ~ 0.24 year). For lack of astronomical data we 
set 3k 2 /(Q = 10" 1 for all other primary-satellite systems considered in Section |ZJ the corresponding 
values of 7 as obtained from (B.3) are given in Table B.2 Of course such values only provide a rough 
guide. 

Table B.2: Values of T (orbital period) and 7 for the systems considered in Section[4j with = 0-1 f° r 

the systems with Jupiter and Saturn as primary. In the third column 7 is computed by using T as time unit, 
whereas the fourth column gives the value of the damping coefficient expressed in years -1 . 





Primary 


Satellite 


T 


T (years) 


7 


7 (years l ) 


E-M 


Earth 


Moon 


2.36x10 s 


7.48X10" 2 


3.75X10" 8 


3.15X10" 6 


S-M 


Sun 


Mercury 


7.60x10 s 


2.41X10- 1 


3.24xl0~ 8 


8.46xl0~ 7 


J-G 


Jupiter 


Ganymede 


6.18x10 s 


1.96xl0~ 2 


1.91xl0- s 


2.48xl0~ 2 


J-I 


Jupiter 


lo 


1.53x10 s 


4.84xl0~ 3 


1.71xl0~ 4 


5.51xl0~ 2 


S-E 


Saturn 


Enceladus 


1.18x10 s 


3.75xl0~ 3 


6.26xl0~ 4 


1.05 


S-D 


Saturn 


Dione 


2.36x10 s 


7.49X10" 3 


1.71xl0~ 4 


1.44X10 -1 
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To obtain the value of e one can use the formula 



h 



2 V p 



3 (Mo\ 
\M )' 



(B.4) 



for the equatorial deformation; see Appendix [C| H ere h% is the tidal Love number (/12 ~ 2ki [HT]). 
while the other constants are as defined after (B.3). If we are interested only in orders of magnitude, 
we can express the equatorial deformation according to (B.4) with hi = 1 for the systems with 
Jupiter or Saturn as primary. Then, by inserting the values of R, p, Mq, M listed in Table B.l| into 



and (B.4) gives 7 



(B.4) and using (B.3) to compute e, we obtain the values in Table 4.1 A comparison between (B.3) 
C e, with C ~ 0.05 (taking the values of the constants ki, £, Q and /12 in the 
0.04 for E-M and C rj 0.055 for S-M). 



literature gives C 

Note that the values so obtained for e are lower than those usually assumed in the literature: 
compare for instance the values e = 2.3x 10" 4 for E-M [58J and e = 1.5xl0" 4 for S-M [3| fT5\ p^l [Taj [27] 

7 and e 



with the corresponding values e = 6.75 x 10 



8.11 x 10 7 in Table 4.1 However, as discussed 



in Section [5j if one does not insist at looking only at the present structure of the satellite, then all its 
evolution plays a relevant role. So, one has to take into account that in the past, when the satellite 
was more fluid, because of the lower value of viscosity, not only the friction was smaller, but also the 
deformation was bigger and hence the coupling e was larger; see also the comments in Section [6j 



C The equatorial deformation 

Consider a homogeneous celestial body S of mean radius R coated by an ocean of depth h > 0, not 
too small. Let P be the centre of attraction. Denote by M and Mq the respective masses and assume 
that the motion of the two celestial bodies about their centre of mass be circular uniform. Let p be 
the distance between the two celestial bodies S and P, with p 3> R 3> h. Assume, for simplicity, that 
S rotates about an axis orthogonal to the plane of the orbit and that the ocean density is negligible 
with respect to the core assumed to be rigid. The discussion below is essentially taken from [36J. 

The distance pc of the centre of mass C from the centre of S is such that pq(Mq + M) = M^p. 
Moreover, if lot denotes the angular velocity of revolution of the two celestial bodies and k is the 
universal gravitational constant, one has = k(M + Mo) by Kepler's third law. Let n be a 

unit vector out of the surface of S and note that, imagining the observer standing on the frame of 
reference rotating around C with angular velocity u, so that the axis from P to S has a fixed unit 
vector i, the potential (gravitational plus centrifugal) energy in the point along the direction n at 
distance r from the centre of S has density V = V5 + Vp + V c { , where 

V s = -k— , V p = -k ^° — , Vc ,= l -u 2 {p 2 c + r 2 -2pcrco^), (C.l) 

r [p^ + r z — zprcosip) 1 / 2 2 

if costp := i ■ n. Expanding V in powers of r/p one finds 

M fr\ 2 f3M 2 , ( 1 



V P + V c{ = -k— I -J UliZ ^ + 2J +const '' (°' 2 ) 

because the linear terms cancel out in virtue of Kepler's third law. Therefore the equation of the 
equipotential surface is 



r \p J V 2 M 2 
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If one writes r = R(l + 8(ip)), with 5(ip) = 5q + 5P 2 (cos where P 2 (z) = (3z 2 — l)/2 is the second 
Legendre polynomial [25J and 5$ is such that the volume of the body P is the same as the volume of 



a sphere of radius R, then (B.3) gives 



S(i/}) = 5P 2 (cos^), 



0. 



Mo 
M ' 



(C.4) 



If the core of the satellite is rigid but the ocean density a is not negligible, e.g. it is equal to 
the core density cr c , then one has to take into account that the tide will modify the potential Vs at 
the site of coordinates r, ip. We make the Ansatz that the equipotential surface is still described as 
r = P(l + 5P2(cos ?/>)), possibly for a different constant 5. Then the density Vs will be 



3kM 
4ttP 3 



7T /'27Z 

sin a da / dip 
n Jo 



R(l+5P 2 (cosa)) 



p 2 dp 



(r 2 + p 2 — 2rp(cos ip cos a + sin ip sin a cos V 7 )) 1 ^ 2 



(C.5) 



By expanding the integrand into Legendre polynomials and using the orthogonality properties of the 
polynomials one finds (see [36], §4.3 for details) 



-kM 



V s 



2R 3 
1 3R 2 



+ g ^3 SP2 (cos ip) 



-kM [ - + -^5P 2 (cosip) 
r 5 H 



r < R, 
r > R. 



(C.6) 



By expanding (C.6) in powers of r/p and summing the leading orders to (C.2), one finds that the 

(C.7) 



equation of the equipotential surface becomes 
kM ( p I \ 2 



P 



R 



1- -<5P 2 (cosV) + 



M 



M 



P 2 (COS'0) 



const., 



up to higher order corrections in R/p. Hence if we look for the constant potential surface we find 

3 



<5(V) = <5P 2 (cosV>), 



. 5/PVMo 
2\p) ~M- 



(C.8) 



which replaces the previous (C.4). The tidal deformation at the surface of the ocean, using the 



notations common in celestial mechanics, can be written as /i2C-fMcos ip), so that the maximal tidal 
excursion is 



h 



-M, C = R 



Mo 
M : 



h 2 



(C.9) 



The number /i 2 is called the tidal Love number. More generally h 2 depends on the detailed structure 
of the satellite, so far supposed to have uniform density cr = a c . If on the contrary cr 7^ cr c , then, 
denoting by r c the shape of the core boundary (while r is the shape of the external ocean surface) , 
one can make once more the Ansatz that the deformations be such that r = R(l + 5P 2 (cos , 0)) and 
r c = i? c (l + <5'P 2 (cos , 0)), where P c is the mean radius of the core and 5,5' are two constants to 
be determined by imposing that the ocean surface is equipotential and balancing the forces acting 
on the core boundary. The latter can be performed by considering the pressures acting on the core 
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boundary due to the elastic forces within the core and the loaded terms caused by the ocean and 
core tide. This leads to two relations involving 5,5' (see [36], §4.4) 



Cc_ 



2oo 
5 a c 



+ 



R c 



1 



Go 
(Tr 



07 



5 Cc 
2 R r 



5' 

3 <T 



R 



1 



0V 



(C.lOa) 



(C.lOb) 



where = (M/M c )<^ (with M c being the mass of the core) and the effective rigidity \x is a dimension- 
less quantity proportional to the rigidity of the core. For instance, if a Q <C cr c , we can approximate 



Rc5' 



5 Cc 



21 + /2 



iM « Cc 



+ 



3 /it 



1 



R J l + jl 



In particular in the limit of high rigidity (// 3> 1) then i? c (5' ~ and R5 ~ (i?/i? c ) 4 Cc = C ( m 
agreement with (B.4)), so that the core deformation becomes very small, that is the core is essentially 
undeformed. 



D Threshold values for the spin-orbit model 



We reason as done for the driven cubic oscillator in [5]- We consider (4.1) with 7 = eC and write it 
as the first order differential equation 

l' = y ' (D.l) 

\y = -eG(9,t)-eC(y-l). 

Then we look for a solution z(t) = (0(t), y(t)) in the form of a power series in e, that is z(t) = z^°\t) + 
ez^\t)+e 2 z^\t) + where z®\t) = (9 + ut,u), with oj = p/q, and z^{t) = {0( k \t), yW(t)) to 
be determined by imposing that z{t) be periodic in t with period 2nq. 
A first order analysis gives 



(0(D= y(D : 

\yW = -G(0 Q + wt,t)-C(u-l). 
By introducing the Wronskian matrix 



(D.2) 



Hi/; (J (D-3) 



we can write z^ (t) as 

C™ S) - H,(f) Q + w ' w i* dT w U* + °> - c ( U -d) • < d ' 4 » 

with to be fixed. Then we obtain 

0« (t) = 0« + y (1 W - / dr /" T dr' [G(0 O + wr', r') + C (w - 1)] , (D.5) 

Jo Jo 
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whereas y^(t) = 9^\t). For (D.4) to be periodic we have to require first of all that 

M(0 O ) 



I r^nq 



/ dt[G(9 + ut,t) + C(u)-l)] = 0, 
Jo 



(D.6) 



then fixi/W in such a way that 



2vrg 



I q dt [ dr [G(6 + ujt, r) + C (w - 1)] 



(D.7) 



while will be fixed to second order by requiring that also 9^ 2 \t) be periodic. 



Using that G(9,t) = dgg(9,t), with g(9,t) given by (4.3), and inserting (D.6) into (D.5) leads to 



1 x P 7 ™? 

)2a k / dt sin(20 o + 2ut - kt) = C (u - 1) 

2nq t£ Jo 



and hence 



2flfe( p ) sin 20 o = - l), fe(p) 



2p 
3 ' 



(D.8) 



(D.9) 



Since / only for = — 3, . . . , 7, k ^ 0, as (4.2) implies, w = p/g is either integer (and 
w G {-1,1,2,3}) or half-integer (and u G {-3/2,-1/2,1/2,3/2,5/2,7/2}). If we confine ourselves 
to positive u, we see that (D.9) fixes 0q provided 



\C\ < C (p/q) : 



_ 2a fc ( p )g 



IP' 



(D.10) 



In particular (D.10) is always satisfied for oj = 1, so that the 1:1 resonance always exists, while for 
the other values of oj we obtain (4.4). 



E Time evolution of friction for the spin-orbit model 

The mechanism of capture into resonance has been studied by several authors starting with the 
theory of capture into the 3:2 resonance of Mercury [231 El E7]. Usually the friction is considered 
either periodic or just not depending on time. Here we regard the friction as not periodic in time 



and given by (2.2), that is starting from a initial very small value, then slowly increasing in order 
of magnitude, until the satellite has completely solidified: such a situation seems possible in the 
formation of a satellite or planet. At the beginning, the satellite can be considered in a fluid state; 
however the dissipation due to tidal effects becomes more and more sensible due to the cooling and 
the resulting increase of viscosity and, eventually, it settles at the final present value: the time over 
which the entire process takes place is called the solidification time and will be denoted by T$. 

We stress that, in the model we are considering, we assume that the satellite has first stabilised 
in its orbit around the primary and then modifies its spinning velocity. Of course the exact evolution 
of the satellite dynamics is still debated and no theory is universally agreed upon. In what follows, 
we shall ignore the model-dependent details of such an evolution. So, for instance, if we accept that 
in some stage of the history of Mercury large quantities of its mantle material have been removed 
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|60| IB], for the purposes of our argument it is not important whether this occurred before Mercury 
attained its final orbital motion or after that event. 

Suppose that friction is essentially due to tidal effects on an originally entirely fluid fast rotating 
satellite (that is with rotation frequency 9 at least a few times larger than the present-day orbital 
frequency ojt) evolving toward a solid body. Assuming that the dynamics is described by (4.1) since 
the beginning could appear contradictory with the fact that originally the satellite was essentially 
fluid, since (4.1) deals with a rigid body with given moments of inertia. On the other hand, as we 
shall argue, friction becomes really effective only when the viscosity has attained high values and, 
when this occurs, the shape of the satellite can be considered close to its final state. One could 
also imagine to study the deformations of the satellite during its evolution, but of course this would 
make the analysis much more complicated; see for instance [J], where asymptotic stability of the 1:1 
resonance is obtained for a deformable body with high rigidity (see also [50]). In other words, using 
the spin-orbit model is justified except possibly for the very early stage of the satellite evolution 
(which are not interested in). 

In the early history of the solar system one can assume that the viscosity is rather low, not too 
different from that of the water, which equals 10 -2 poises (CGS units). The solidification process 
is very complicated, and although it has been extensively studied in the literature, especially in the 
case of the Earth and the Moon (for the obvious reason that it is much easier to compare the results 
obtained by theoretical methods and numerical simulations with the experimental data), still there 
are many unsolved issues. See for instance [51 ] for a review of recent results on the Moon. Usually 
one assumes that originally all satellites and planets were completely or almost completely molten, 
according to the so-called magma ocean hypothesis; in the case of the Moon see [57] and especially 
|51j and the references quoted therein (see also [55] in the case of Earth) . Most of the celestial body 
crystallised, from the bottom to the top, following a sequence determined by chemical composition 
and pressure [52], [291 EB] leaving only a layer of very fluid magma close to the surface of the satellite. 
The outer liquid part eventually disappeared, through a further solidification process from the core- 
mantle boundary to the surface, and it is irrelevant in any case to the damping, because of its very 
low viscosity. 

Timing of formation of satellites and for the cooling of the magma is mostly deduced from the 
study of rock samples: in the case of the Moon of course this is much easier |51 | llO | [T7 1 l38|. while the 
results are far from being conclusive in the case of Mercury (however see, for instance, [531 111] ). In 
any case, it is not unlikely that the solidification time of Mercury is of the same order of magnitude 
of that of the Moon; we can also mention that the theory has been proposed that the two celestial 
bodies may have had similar origin |59| 160] . There is strong evidence that the solidification time T$ 
for the Moon is of order of 10 s years [101 138] . Thus, for the reasons given above, we take this as the 
order of magnitude of the solidification time of all the satellites considered in Section [4] 

For most of the solidification process the ocean magma is maintained above its solidus [HUB]. So 
it is reasonable that the viscosity has increased as an effect of the cooling of the satellite: eventually 
most of the satellite is almost solidified and its viscosity is very large, say of the order of that of the 
mantle, which is almost solid (hence about 10 24 poises, much higher than the viscosity of the terrestrial 
magma which ranges from 10 4 and 10 12 poises \42\ [30]). Of course the viscosity strongly depends 
on temperature, which in turn decreases in time during the process of cooling of the satellite. One 
could look at the literature for profiles of temperatures versus time |55| |2T] or for the dependence 
of viscosity on temperature in fluids [49J. A detailed discussion on viscosity evolution during the 
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solidification of the satellite stands as a very hard problem, and in fact such an issue is not widely 
studied [13 [181 [51] ; see however [46J, where the despinning of Saturn's satellite Iapetus is studied 
and an Arrhenius law is assumed to describe the temperature dependence of the viscosity, and |47j , 
where the despinning of Mercury is related to the thermal evolution of the planet. 

As far as the satellite can be considered essentially fluid, the damping coefficient 7 has to be 
formed with the following physical quantities: the viscosity 77 of the magmatic fluid constituting the 
satellite, the mass M of the satellite, the equatorial deformation 2h due to the tides and the angular 
velocity ujt- Starting from the rescaled equation (4.1), one expects 7 oc rjR 2 h/I z , with I z oc MR 2 ; 
then, by rescaling time t — > t/uiT, an appropriate choice for 7 turns out to be 7 ~ -qh/Mux- 
The evolution of magma from the initial melt passes through the formation of cumulate rocks and 
fractional crystallisation, leading to the reduction of the outer liquid layer with low viscosity and the 
accretion of the internal, highly viscous core; for a more detailed discussion see for instance [51J in the 
case of the Moon. Therefore, when the solidification process attained a high stage of advancement, a 
different model for the friction has to be taken. If the satellite is essentially solid, with a thin external 
fluid layer (ocean), one can neglect the fluid part, because of its very low viscosity, and concentrate 
on the solid core. One can use the analysis in Appendix [Bj in the limit case a Q <C a c and very high 
rigidity fi, so that the core deformation in (C.lOb) is very small. Also in this case, one can express 
7 in terms of the involved physical quantities and, again on the basis of dimensional arguments, set 
7 ~ /j/i/Mwy. Therefore, to sum up, friction increases with viscosity up to a certain value. Once 
such a value is reached, at which the satellite can be considered essentially a solid with very high 
rigidity. 

The despinning time Tjj represents the time which the satellite needs to be really attracted into 



a resonance. Hence (see Appendix |A|) Td = 0(7 _1 ), where 7 is the quantity appearing in (4.1). Of 
course, we would want that the solidification time be larger or at least comparable with despinning 
time. This is the case if one can assume T5 ~ 10 8 years and Td ~ I/7, with 7 as in Table B.2 (which 
yields Td « 10 6 years). 



F Numerical details 



The numerical results have been found by running a variety of computer programs which implement 
different algorithms. The main algorithms used were (i) a standard Runge-Kutta integrator with 
automatic step-size control [44]; (ii) the Bulirsch-Stoer algorithm [44], which extrapolates the step- 
size to zero and is suitable for high-accuracy computation; and (hi) a numerical implementation of 
the Frobenius method. Of these, (ii) is the slowest, but serves to confirm the results obtained from 
the other methods. Most of the data in this paper came from (i) and (iii), of which (iii) is the fastest. 
This is because the use of series often enables large time steps to be taken. On the assumption that 
the solution to the differential equation around a point t = to can be expressed as a power series 
in t, we obtain a (nonlinear) recursion formula for the coefficients in this series. In principle, as 
many terms as desired can be computed, and in practice about 25 worked well. For a given initial 
condition, this series enables us to compute x and x for \t — to\ < R, for some R that depends on the 
desired accuracy, the initial conditions, and to- The size of the step, t — to, is chosen as the one that 
makes the absolute value of the right-hand side of the differential equation, which should of course 
be zero, smaller than some tolerance: 10 -12 was used here. Typical step sizes ranged between about 
0.29 and 10.0, with 0.7 being chosen about 50% of the time. Even the smallest step size is many 
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times larger than that which would be used in a Runge-Kutta implementation. 

The initial conditions are chosen randomly and are uniformly distributed inside the square Q. 
Since we seek detailed estimates of the relative areas of the basins of attraction, at least up to the 
first decimal digit whenever possible, we have to take many initial conditions: certainly 1000 initial 
conditions, as in [T3] is not enough. For the data to be reliable to the first decimal place, with a 
95% confidence level, we found that 1 000 000 initial conditions need to be taken inside the square Q 
defined in Section[2} A statistically rigorous justification for this confidence level, given the number of 
initial conditions, can be found in [56, chapter 9]. Thus, for the relative areas in Table 2.3, we used up 
to 1 000 000 initial conditions except for smallest values of 7 (500 000 initial condition for 7 = 0.00005, 
150 000 initial condition for 7 = 0.00001 and 50 000 initial conditions for 7 = 0.000005). Analogously 
we considered 1000 000 initial conditions for e = 0.5 (Table pTl]) and 500 000 initial conditions for 



e = 0.01 (Table 5.2). Also in Section [3j we considered as many initial data as possible: 1000 000 
initial conditions for larger 7 (7 = 0.015), 500 000 for 7 = 0.005 and 250 000 for the smallest value 
7 = 0.0005. As a general rule, the smaller 7 is, the longer the integration time and hence, by taking 7 
smaller, we also need to diminish the number of initial conditions, in order that the computation time 
does not become prohibitively long. However, the error on the relative areas becomes larger. That 
said, the general scenario described in Section [3] seems clear and well supported by the numerics. 

As already noted in [JH], for 7 very small, the basins of attraction are spread out over the entire 
phase space and become very sparse. Thus, if one wants to detect which basin a given initial condition 
belongs to, very high numerical precision is needed. 

Finally, the integration time T[ nt must be chosen in such a way that all trajectories reach the 
attractor (within a reasonably fixed accuracy). For instance, one can take Ti n t = iV/70, with iV = 20. 
Therefore, in order to investigate the dynamics for very small values of the damping coefficient, T; n t 
has to be very large. 

The conclusion is that we have to follow the trajectories of a large number of initial conditions, 
for very long times and with very high accuracy. Of course, this is at odds with obtaining results 
within a a reasonable time, so that we need to reach a compromise. This has led us to the choice 
described above. 



References 

[1] Y. Abe, Thermal and chemical evolution of the terrestrial magma ocean, Phys. Earth Planet. Int. 100 
(1997), 27-39. 

[2] E.Yu. Aleshkina, Synchronous spin-orbital resonance locking of large planetary satellites, Solar Syst. Res. 
43 (2009), no. 1, 71-78. 

[3] J.D. Anderson, G. Colombo, P.B. Esposito, E.L. Lau, G.B. Trager, The mass, gravity field, and ephemeris 
of Mercury, Icarus 71 (1987), 337-349. 

[4] D. Bambusi, E. Haus, Asymptotic stability of spin orbit resonance for the dynamics of a viscoelastic 
satellite, Preprint, 2010. 

[5] M.V. Bartuccelli, A. Berretti, J.H.B. Deane, G. Gentile, S. Gourley, Periodic orbits and scaling laws for 
a driven damped quartic oscillator, Nonlinear Anal. Real World Appl. 9 (2008), no. 5, 1966-1988. 



27 



[6] M.V. Bartuccclli, J.H.B. Deane, G. Gentile, S. Gourley, Global attraction to the origin in a parametrically 
driven nonlinear oscillator, Appl. Math. Comput. 153 (2004), no. 1, 111. 

[7] M.V. Bartuccelli, G. Gentile, K.V. Georgiou, On the dynamics of a vertically driven damped planar 
pendulum, R. Soc. Lond. Proc. Ser. A Math. Phys. Eng. Sci. 457 (2001), no. 2016, 3007-3022. 

[8] W. Benz, A. Anic, J. Horner, J.A. Whitby, The origin of Mercury, Space Sci. Rev. 132 (2007), 189-202. 

[9] G.D. Birkhoff, Dynamical systems, American Mathematical Society Colloquium Publications Vol. IX, 
Providence, R.I., 1966. 

[10] A. Brandon, A younger moon, Nature 450 (2007), 1169-1170. 

[11] S.M. Brown, L.T. Elkins-Tanton, Compositions of Mercury's earliest crust from magma ocean models, 
Earth Planet. Sci. Lett. 286 (2009), 446-455. 

[12] T.J. Burns, C.K.R.T. Jones, A mechanism for capture into resonances, Phys. D 69 (1993), no. 1-2, 
85-106. 

[13] A. Celletti, L. Chierchia, Measures of basins of attraction in spin-orbit dynamics, Celestial Mech. Dynam. 
Astronom. 101 (2008), no. 1-2, 159-170. 

[14] Sh.N. Chow, J.K. Hale, Methods of bifurcation theory, Grundlehren der Mathematischcn Wissenschaften 
Vol. 251, Springer, New York-Berlin, 1982. 

[15] A. CM. Correia, J. Laskar, Mercury's capture into the 3/2 spin-orbit resonance as a result of its chaotic 
dynamics Nature 429 (2004), 848-850. 

[16] A. CM. Correia, J. Laskar, Mercury's capture into the 3/2 spin-orbit resonance including the effect of 
core-mantle friction, Icarus 201 (2009), 1-11. 

[17] V. Debaille, A. D. Brandon, Q. Z. Yin, B. Jacobsen, Coupled 142 Nd- 143 Nd evidence for a protracted 
magma ocean in Mars, Nature 450 (2007), 525-528. 

[18] L.T. Elkins-Tanton, Linked magma ocean solidification and atmospheric growth for Earth and Mars, 
Earth Planet. Sci. Lett. 271 (2008), 181-191. 

[19] U. Feudel,C. Grebogi, Why are chaotic attractors rare in multistable systems, Phys. Rev. E 91 (2003), 
134102, 4 pp. 

[20] U. Feudel,C. Grebogi, B.R. Hunt, J.A. Yorke, Map with more than 100 coexisting low-period periodic 
attractors, Phys. Rev E 54 (1996), 71-81. 

[21] P.E. Fricker, R.T. Reynolds, A.L. Summers, Possible thermal history of the Moon, in The Moon, Eds. 
H.C Urey and S.K. Runcorn, Reidel, Dordrecht, pp. 384-391, 1972. 

[22] G. Gentile, M.V. Bartuccclli, J.H.B. Deane, Bifurcation curves of subharmonic solutions and Melnikov 
theory under degeneracies, Rev. Math. Phys. 19 (2007), no. 3, 307-348. 

[23] P. Goldreich, S. Peale, Spin-orbit coupling in the solar system, Astronom. J. 71 (1966), no. 6, 425-438. 

[24] P. Goldreich, S. Sotcr, Q in the Solar system, Icarus 5 (1966), 375-389. 

[25] I.S. Gradshteyn, I.M. Ryzhik, Table of integrals, series, and products, Sixth edition, Academic Press, 
Inc., San Diego, 2000. 



28 



[26] J. Guckenheimcr, Ph. Holmes, Nonlinear oscillations, dynamical systems, and bifurcations of vector fields, 
Applied Mathematical Sciences Vol. 42, Springer, New York, 1990. 

[27] T. Van Hoolst, F. Sohl, I. Holin, O. Verhoeven, V. Dchant, T. Spohn, Mercury's interior structure, 
rotation, and tides, Space Sci. Rev. 132 (2007), 203-227. 

[28] A. Khan, K. Mosegaard, J. G. Williams, P. Lognonnc, Does the Moon possess a molten core? Probing the 
deep lunar interior using results from LLR and lunar prospector, J. Geophys. Res. 109 (2004), E09007, 
25 pp. 

[29] Th. Kleinc, H. Palme, K. Mezger, A.N. Halliday, Hf- W chronometry of lunar metals and the age and 
early differentiation of the Moon, Science 310 (2005), no. 5754, 1671-1674. 

[30] D.L. Kohlstcdt, M.E. Zimmerman, Rheology of partially molten mantle rocks, Ann. Rev. Earth Planet. 
Sci. 24, (1996), 41-62. 

[31] N.N. Krasovsky, Stability of motion. Applications of Lyapunov's second method to differential systems 
and equations with delay, Stanford University Press, Stanford, Calif. 1963. 

[32] J. Laskar, O. Neron de Surgy, On the long term evolution of the spin of the Earth, Astronom. Astrophys. 
318 (1997), 975-989. 

[33] A.J. Lichtenberg, M.A. Lieberman, Regular and chaotic dynamics, Applied Mathematical Sciences Vol. 
38, Springer, New York, 1992. 

[34] G.J.F. MacDonald, Tidal friction, Rev. Geophys. 2 (1964), no. 3, 467-541. 

[35] W. Magnus, S. Winkler, Hill's equation, Dover Publications, Inc., New York, 1979. 

[36] CD. Murray, S.F. Dermott, Solar system dynamics, Cambridge University Press, Cambridge, 2001. 

[37] A.I. Neishtadt, On probabilistic phenomena in perturbed systems, Selecta Math. Sov. 12 (1993), no. 3, 
195-210. 

[38] A. Nemchin, N. Timms, R. Pidgeon, T. Geisler, S. Reddy, C. Meyer, Timing of crystallisation of the 
lunar magma ocean constrained by the oldest zircon, Nature Geosci. 2 (2009), 133-136. 

[39] J. Palis, A global view of Dynamics and a conjecture on the denseness of finitude of attr actors, Asterisque 
261 (1995) 335-347. 

[40] J. Palis, A global perspective for non-conservative dynamics, Ann. Inst. H. Poincare Anal. Non Lineairc 
22 (2005), no. 4, 485-507. 

[41] S.J. Pcalc, The free precession and libration of Mercury, Icarus 178 (2005), 4-18. 

[42] C.C. Plummer, D. McGeary, D.H. Carlson, Physical geology, McGraw-Hill, 2003. 

[43] M. Porter, P. Cvitanovic, A perturbative analysis of modulated amplitude waves in Bose-Einstein con- 
densates, Chaos 14 (2004), no. 3, 739-755. 

[44] W.H. Press, S.A. Tcukolsky, W.T. Vetterling and B.P. Flannery, Numerical Recipes in C, Cambridge 
University Press, Cambridge, UK, 1992. 

[45] M.E. Pritchard, D.J. Stevenson, Thermal implications of a lunar origin by giant impact, in Origin of the 
Earth and Moon, Eds. R.M. Canup and K. Righter, University of Arizona Press, pp. 179-196, 2000. 



29 



[46] G. Robuchon, G. Choblet, G. Tobie, O. Cadck, C. Sotin, O. Grasset, Coupling of thermal evolution and 
despinning of early Iapetus, Icarus 207 (2010), 959-971. 

[47] G. Robuchon, G. Tobie, G. Choblet, O. Cadek, A. Mocquet, Thermal evolution of Mercury: implication 
for despinning and contraction, 40th Lunar and Planetary Science Conference (2009, The Woodlands, 
Texas), id. 1866. 

[48] Ch.S. Rodrigues, A.P.S. dc Moura, C. Grebogi, Emerging attractors and the transition from dissipative 
to conservative dynamics, Phys. Rev. E 80 (2009), no. 2, 026205, 8 pp. 

[49] Ch.J. Seeton, Viscosity-temperature correlation for liquids, Tribology Lett. 22 (2006), no. 1, 67-78. 

[50] A.V. Shatina, Evolution of the motion of a viscoelastic sphere in a central Newtonian field, Cosmic Res. 
39 (2001), no. 3, 282-294. 

[51] C.K. Shearer et at, Thermal and Magmatic Evolution of the Moon, Rev. Mineral. Geochem. 60 (2006), 
365-518. 

[52] G.A. Snyder, L.A. Taylor, C.R. Neal, A chemical model for generating the sources of mare basalts: 
combined equilibrium and fractional crystalization of the lunar magmasphere, Geochim. Cosmochim. Acta 
56 (1992), 3809-3823. 

[53] S.C. Solomon et at, The MESSENGER mission to Mercury: scientific objectives and implementation, 
Planet. Space Sci. 49 (2001) 1445-1465. 

[54] T. Spohn, F. Sohl, K. Wieczerkowski, V.Conzelmann, The interior structure of Mercury: what we know, 
what we expect from Bepi-Colombo, Planet. Space Sci. 49 (2001) 1461-1470. 

[55] D.L. Turcotte, On the thermal evolution of the Earth, Earth Planet. Sci. Lett. 48 (1980), 53-58. 

[56] R.E. Walpole, R.H. Myers and S.L. Myers, Probability and Statistics for Engineers and Scientists, Prentice 
Hall, New Jersey, US, 1998. 

[57] P. Warren, The lunar magma ocean concept and lunar evolution, Ann. Rev. Earth. Planet. Sci. 13, (1985), 
201-240. 

[58] J.G. Williams, X.X. Ncwhall, J. O. Dickey, Lunar moments, tides, orientation, and coordinate frames, 
Plan. Space Sci. 44 (1996), 1077-1080. 

[59] M.M. Woolfson, The origin and evolution of solar system, Graduate series in Astronomy, Institute of 
Physics Publishing, Bristol, 2000. 

[60] M.M. Woolfson, On the origin of planets, Imperial College Press, London, 2011. 

[61] C.Z. Zhang, Love numbers of the Moon and of the terrestrial planets, Earth, Moon and Plan. 56 (1992), 
193-207. 



30 



